#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <math.h>

void seqScan(double *X, unsigned int n, unsigned int offset, unsigned int stride);

void seqAdd(double *X, double val, unsigned int size, unsigned int offset, unsigned int stride);

// chunkBits is the number of bits for a block, presumable fits in cache.
void hybridScan2(double *X, int n, int chunkBits) {
	int i, j;
	unsigned int d, logn, numElements;
	unsigned int numEl, block, blockBits;
	unsigned int stride = 1;

	logn = (int)ceil(log2(n));

	// shift the 1 into the chunkBits+1 position	
	// always deal with these many in each iteration
	numElements = 1 << chunkBits;  
	
	for (d = 0; d < logn; d += chunkBits) {
		  // perform seq scan on 2^chunkBits of entries at a time
		blockBits = d + chunkBits;

		  // limit the block size
		blockBits = (blockBits > logn ? logn : blockBits);   
		block = 1 << blockBits;

		// do blockwise sequential scan
#pragma omp parallel for shared(X, numElements, n, stride, block) \
	private(numEl) schedule(guided)
		for (i = 0; i < n; i += block) {
			numEl = (i + block) > n ? (n - i) / stride : numElements;
			seqScan(X, numEl, i + stride - 1, stride);
		}

		stride <<= chunkBits;
	}
	
	int leftOvers = logn % chunkBits;
	// finally do blockwise sequential update at the very end.
	block = 1 << (logn - leftOvers);
	for (d = logn - leftOvers; d >= chunkBits; d -= chunkBits) {
		// perform update 
		stride = block >> chunkBits;
#pragma omp parallel for shared(X, numElements, n, block) \
	private(numEl) schedule(guided)
		for (i = block; i < n; i += block) {
			// need to replae >= with >
			numEl = (i + block) > n ? (n - i) / stride : numElements - 1;
			seqAdd(X, X[i - 1], numEl, i, stride); // and use the previous value
		}
		block >>= chunkBits;
	}
}
